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Abstract 

Self-organization on crystal surface is studied as a two dimensional spinodal 
decomposition in presence of a surface stress. The elastic Green function is 
calculated for a (001) cubic crystal surface taking into account the crystal 
anisotropy. 

Numerical calculations show that the phase separation is driven by the in- 
terplay between domain boundary energy and long range elastic interactions. 
At late stage of the phase separation process, a steady state appears with dif- 
ferent nanometric patterns according to the surface coverage and the crystal 
elastic constants. 
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I. INTRODUCTION 



Self-organization (SO) on solid surface is an efficient mean for growing nanostructures 
with regular sizes and spacings. The models proposed by Marchenko [0] and Vanderbilt et 
al. |2||| are the basis of the theoretical framework to understand the SO phenomenon. They 
enhanced the interplay between the long range elastic interaction yielded by the underlying 
crystal and the domain boundary energy Indeed, the former is minimum when two surface 
defaults are separated by a distance as large as possible while the latter is minimum when 
only one compact domain appears onto the surface. So when these two ingredients are 
present (see below for experimental descriptions), the surface ground state structure should 
balance the aforementioned interactions. The purpose of was to evaluate this structure 
in different cases. 

Assuming that two phases, called A and B, coexisting onto a crystal surface, have differ- 
ent intrinsic stresses, the authors of Refs. showed that a state which consists of stripes 
domains occupied alternatively by A and B lowers the energy with a period selection. This 
period depends on the crystal stress energy compared with the domain boundary energy 
and it varies exponentially with the ratio of those two quantities. In order to perform an 
analytical surface elastic Green function calculation, Marchenko and Vanderbilt both as- 
sumed either a crystal anisotropy along one direction of the surface or the anisotropy of the 
intrinsic stresses. Those works specially addressed the cases of corrugated crystal surfaces 
and crystal surface reconstruction with broken symmetry (see Ref. for recent relevant 
experimental analysis). 

Comparing the stability of different periodic domains in an isotropic 2-dimensional dipo- 
lar model, Vanderbilt et al. @ proved that at intermediate coverage, 9$ > 0.28, the stripe 
structure is the optimal candidate for the 2-dimensional system ground state while at low 
coverage, 9q < 0.28, the droplets structure is more stable. As the dipole-dipole interaction is 
similar to the elastic interaction, i.e., it decreases as the inverse of the distance to the power 
three, these Vanderbilt's results hold for the solid surface SO. It shows that the anisotropy of 
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the intrinsic stress surface, assumed in Refs. |],|2| is not essential in the SO process since the 
patterning occures even if the surface is isotropic with no symmetry breaking. The studies 
of Refs. [0-0] emphasize the main physical ingredients which drive SO, i.e. , the interplay 
between domain boundary energy and elastic interaction energy. 

Recent analysis of chemisorbed mono-layers on (001) Copper surfaces, via Scaning Tun- 
neling Microscopy (STM) f)|-|8[ and Spot Profile Analysing Low Energy Electron Diffraction 
(SPA-LEED) |§ show mesoscopic morphologies different from droplet or stripe structures. 
For the N/Cu(001) case, the Nitrogen is chemically adsorbed on a (001) Copper crystal sur- 
face, and agregates within square-shaped islands that may arrange either in 1-dimensional 
rafts at low Nitrogen coverage or in a 2-dimensional array at intermediate coverages. The 
experimental works mentioned above |7],[|] motivated the present study. 

Here we propose to describe the SO kinetics on solid surface as a 2-dimensional spinodal 
decomposition. In this standard theory, we include the stress energy due to the underlying 
crystal and therefore the surface elastic Green function is calculated. The 2D Cahn-Hilliard 
equation which drives the surface diffusion, is integrated with computer means. Our method 
allows us to capture the SO kinetics together with the elastic anisotropy due to the crystal 
symmetries. The latter feature is proved to play a role in the nanometric arrangements. For 
simplicity, our study is focused on a (001) cubic crystal surface. 

In the last decade, a similar theory has been developed under the name of "phase field" 
by Khachaturyan et al. for the kinetics of phase transition in alloys [[K]]. The phase field 



name is also used for models of solidification; see |TTj] for recent developments. We choose 
not to use this ambiguous terminology. In the framework of the SO spinodal theory, the 
implied coexistent phases may represent either two types of crystal facet, noted A and B or 
a chemical adsorbed layer A over a clean crystal surface B. In the latter case, we neglect the 
possible layer thickness. 

Our results exhibit a steady state at late stage of the separation process which shows 
different mesoscopic patterns according to the coverage and to the elastic constants of the 
materials. Because of the crystal cubic symmetry, the nanometric morphologies of the final 



state may differ from the one predicted in Refs. jJT^Q] . At low coverage, rafts of either disk 
or square shaped islands appear while at intermediate coverage a branched stripe structure 
occurs. The elastic constants of the material are shown to determine the preferential orien- 
tations of the rafts and of the stripes. If the surface square symmetry is broken, the kinetics 
final state is similar to those predicted in Refs. i.e. , an assembly of regular spaced 

stripes. 

II. DOMAIN BOUNDARY ENERGY 

The covering parameter 9 describes the A-B coexistent phases at the solid surface. These 
phases may differ by their composition or their geometry. Let say that 9 = for B and 
9 = 1 for A. If one assumes that the spatial surface variations of 9 are smooth with respect 
to the atomic scale, a 2-dimensional coarse-graining procedure is thus relevant to represent 
the mesoscopic system state. If d is the size of the elementary coarse-grained surfaces, 
we introduce the mesoscopic quantity 9(r) which is the average, performed over a whole 
elementary surface, centered at position r. The 9(r) variable is the mesoscopic local coverage 
of the surface since it is the A quantity per unit area which is present in the r vicinity. The 
^-average over the entire surface is written 9q and it is assumed to be conserved during the 
system time evolution. 

An inhomogeneous mixing of A and B phases involves an energy increasing because of 
atomic bond breaking at domain boundaries. This energy together with the entropic term 
due to the surface inhomogeneity are capted in the 2-dimentional free Ginzburg-Landau 
(GL) functional energy: 

F -»=^-// s {i0 2+( £ )2|+/w>rfr (21) 

We introduce here the adimensional free energy density / = 16.0 2 (1 — 9) 2 which is a double 
well potential with minima for 9 = and 9 = 1. The hat notation points out the adimen- 
tional quantities. The / form has no direct influence on the mesoscopic structure providing 
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GL functional is invariant with respect to the (001) surface space group. Let note g(0), the 
continuous Hamiltonian including both the 2D gradient term and the free energy density /. 

The Fq and 7 scalars are respectively the free energy density constant and the amplitude 
of the gradient term that both must be adjusted to set the model domain boundary energy 
/ to a realistic value, i.e., around 10 meV/A (see Ref. 0). At phase equilibrium, the I 
quantity is given by / = Fq f£ g(0)dl where the integration is performed along a line path 
that goes from inside a B phase domain to inside an A domain [fT^] . It is easy to see that 
the previous path integral is overvalued by the product of the A-B interface width times 
the / maximum value, i.e., f{9 = 0.5) = 1. Fixing 7 = 20 d 2 , which insure the 9{r) 
space variation smoothness, the A-B interface width is then around 5 . d, at equilibrium. 
Therefore, to obtain / m 10 meV/A, F .d must be fixed to 2 meV/A (or 3.2 10~ 12 J/m). 
As our investigations are focused over the nanometer scale, we choose to set d = 1 nm and 
thus Fq = 3.2 mJ/m 2 . To study larger space scales as it would be suitable for Silicon which 
may exhibit a 100 nm vicinal period [[|, then it is sufficient to increase the d parameter and 
to adjust subsequently F , eventually changing the domain boundary energy I if necessary. 

For simplicity, we choose to neglect the possible 7 variations with respect to the crystal 
surface direction. Such a feature may be yielded from a crystal step anisotropy but here our 
study is focused on the elastic anisotropy due to the crystal symmetries. 

III. SURFACE GREEN FUNCTION 

An important stage of this work consists in the calculation of the cubic crystal surface 
elastic Green function. Let note P(r) a surface external force at position r. In cartesian 
coordinates, the (001) surface is defined by (£3 = 0) and r = (07, x 2 )- The semi-infinite 
crystal occupies the half space 23 > 0. The surface normal is the n unit vector. The 
mechanic equilibrium condition at the surface is given by: 

(7ij(r,x 3 = 0).nj = Pj(r) (3.1) 
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where rij is a n component and the summation over subscript j is implicit. The crystal 
bulk stress, aij(r,x 3 ) is due to the crystal displacements and these quantities are related 
to each other by the Hooke law: a^j = Xij,k,iduk/dxi. The forth order tensor \ij,k,i gives 
the crystal elastic constants and for a cubic crystal symmetry, this tensor is composed with 



three non zero coefficients [FJJ], namely A 



Cn, X 



' "1,1,3,3 



C12 and A 



1,3,1,3 



A 



i,3,3,i 



C44. 



In case of copper, those coefficients are : Cn = 1.683 I0 11 J/m 3 , C12 = 1. 221 I0 11 J/m 3 and 
C 44 = 0.757 I0 11 J/m 3 (see Physics Handbooks). 

At mechanic equilibrium, the bulk displacements fulfill the Lame equation: 

d 2 u k 



A 



i,j,k,l 



dxjdxi 







(3.2) 



Except constants, the displacement functions are fully determined by the set of equations 
( j3TT| , |3.2|) . As proposed in Ref. [|15[ for the isotropic case, these equations may be solved by 
writing the displacements as 2-dimensional Fourier transforms: 



Ui(r,x 3 ) = J J exp (iQ.r).Ui(Q,x 3 )dQ 



(3.3) 



Note that the Fourier components Ui(Q,x 3 ) depend on both the wave vector Q = (91,52) 
and x 3 . Once the set of Eqs. ( |3.I| , |3.2p is expressed using Eq.(|3.3|), the subsequent differential 
equations involving Ui and their derivatives with respect to x 3 turn out to be linear, so that 
£ti(Q 7^ 0,0:3) have an exponential dependence on the bulk penetration length x 3 : 



Uj(Q ^ 0,x 3 ) = ^/3 jV exp {-ai.x 3 ) 
1 



(3.4) 



where both [3jj and a\ depend on the wave vector Q. As the whole sample is at rest, the P 
average is assumed to be zero which implies Uj(0, x 3 ) = 0. The Eq. ( |3.2| ) is then reduced to 
a linear equation M(Q, ai)(j3jj) = where M(Q, cti) is the 3x3 matrix: 



M(Q,< 



^ C n qf + C u q\ - C AA af (C 12 + C AA )q x q 2 (C 12 + C AA )iq x ai 

(C12 + C 44 )gig2 C\ X q\ + C M q\ - C u af (C 12 + C iA )iq 2 ai 
1 -(C12 + C AA )%q x (Xi -(C12 + C u )iq 2 ai C n af - C u (qf + q\ 



\ 



(3.5) 
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As usual, the non-trivial solutions are such as det (M(Q,«/)) = 0. This yields a third 
degree polynomial equation for the atf parameters that we solve numerically. The only three 
a\ values with positive real parts are physically acceptable. At least one of the on is real, 
the two other roots may be either real or complex conjugated depending on the sign of the 
X parameter: 

X = C n - C12 - 2C 44 (3.6) 

This combination of the elastic constants is related to the elastic anisotropy of the cubic 
crystal (see Ref. |14[]). When x — 0, the crystal is isotropic and the three aj are degenerate. 
This case was addressed in Refs. |14l , |T5|l . One retains as examples that the Copper and Gold 



x's are negative (xcu = —1-0, XAu = —0.5) and the Chromium and Niobium x's are positive 
(XCr = +1-8, XNb = +0.5). 

To each aj corresponds a unique set of three coefficients (Pj t i) that are determined by 
inverting Eq. ( |3.1| ). Noting Pj the Fourier transform of the force component Pj, we write 
the Fourier counterpart of the Eq. (EOT): 



-aiPiA + %A,3 = Pi 

'=(1,2,3) 
'=(1,2,3) 

C 12 (iqiPi,i + iqzPl?) ~ C u aif3i,3 = P3 

'=(1,2,3) 

(3.7) 

Combining the first and second rows with the third rows of the matrix M(Q, ai) (see 
Eq. (|3.5|) ), one easily gets both equations 

A,i = iQiA,a r ',i/"' 

A,2 = «<?2-A,3 ^l,2/0Cl 

(3.8) 

where we write: 



(Cn - C12 ~ CuWi - C 44 (g 2 + g|) 
(Cii - Ci 2 - C 44 )g 2 + C 44 (g 2 2 - a?) 
{Cn ~ Cn - - C 4 M + g|) 

(C n - C 12 - C u )ql + CM - a?) 



(3.9) 



provided the denominators of r^i and T i 2 are non zero. If not, which occures for either 
qi = or q 2 = then it is easy to show that ^3 = and the calculation can be performed 
as well. In addition let us note 



1,3 



-C n ai - Ci 2 (g 2 r Zi i + qlTi^)/ai 



(3.10) 



We rewrite Eq. ( |3.7| ) as a linear equation Njjftip = Pj where the 3x3 matrix N is defined 
as follows: 



^ %(1 - r\i) igi(l-r 2 ,i) igi(l-r 3j i)^ 



N(Q) 



V 



(3.11) 



ig 2 (l-r 1)2 ) ig 2 (l-r 2j2 ) ig 2 (l-r 3i2 ) 

Ti,3 r 2 ,3 r 3j3 

Provided gi 7^ g 2 , this matrix ( |3.11 ) can be inverted which gives the Pi^s as linear functions 
of the Pj(r) Fourier transforms : /3/ >3 = NfJPj. One must distinguish the case q% = q 2 for 
which details are not presented but are simple to deal with. Using Eqs. (|3.8| ), one gets the 
whole set of $j's an d thus Uk(Q,X3 = 0) = Gk,j(Q)Pj(Q) where Gk,j is the surface elastic 
Green function that we write as follows: 

1 

g 2iJ (q) = ^e^;/ 
1 

c 3 ,AQ) = E^ 

(3.12) 

The latest function varies as 1/|Q|. The total elastic energy of the system is given at 
mechanical equilibrium according to Eqs. (|3~T| , |372|) by : 



E el = -1/2 / Poinds (3.13) 

Jx 3 =0 

that gives an analytical expression of E e \ in the Fourier space: 

E d = -l/2[ P^G^dQ (3.14) 

The total energy is in fact given by the sum F = F c h em + E e i (Eqs. ( p7T| ) and ( |3.2| )). 

The P force is induced by the phases intrinsic stress misfit. Let note a A and a B the 
intrinsic stress tensor of the A and B phases, respectively. At position r, the intrinsic stress 
is thus given by a°(r) = a B + 5<j°0(r) where we introduce the tensor 5a° = a A — a B . The 
induced force is simply obtained by deriving a with respect to the surface coordinates, it 
gives: 

Here we propose to distinguish whether there is a strong directional anisotropy such 
as for vicinal surface or no broken surface symmetry. In the former case, if a negligible 
stress variation is induced along a direction which angle with respect to [100] is noted a, 
the P vector must be zero if gradient points this direction. The 5a° tensor coefficients are 
developed to satisfy this condition and one may write P as follows: 

80 80 

Pi = Aj [sin(a]- cos (a) — — ] (3. 16) 

ox % OX2 

where we introduce three constants {Aj}. We restrict our study to a uniaxial stress with 
only one non zero Aj which is A 3 = A. According to our tests, non zero K\ and A 2 involve 
no qualitative change in our results. 

If the surface square symmetry is preserved, the Sa° tensor is invariant under the square 
symmetry group operators. So 5cr° 2 = Sa^i, — ^^32 an d = da^- Again we 

introduce the A parameter setting that Sa^ = = A. The surface is assumed to be a 
perfect plane so the P3 component is neglig ible, i.e., Sa^ = 5a% 2 = 0. Non zero coefficients 
\i = 5a1 2 = £o"2i implies a shear that may be induced by the internal structures of the 
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coexistent phases. To a first step, we choose to ignore such a possibility, i.e., \x << A. Thus 
the force P is simply proportional to the 9 gradient. 

The parameter A is a constant that fixes the amplitude of the surface force and it 
may vary from an experimental case to another. In our model, for each different set of 
crystal elastic coefficients, A is adjusted, i.e., our numerical calculations are performed with 
A = 40 mJ/m 2 (or A = 0.25eV/A 2 ) for the Copper, Gold and Niobium crystal cases and 
A = 57 mJ/m 2 (or A = 0.35eV/A 2 ) for the Chromium case. 



IV. KINETICS OF PHASE SEPARATION 



To modelize the phase separation kinetics, we use a standard spinodal decomposition 



theory which is precisely described in Ref. [12| for the 3-dimensional case. Extension to the 



2-dimensional systems is straighfoward. The time evolution of the conserved quantity 9(r) 
is then driven by the Cahn-Hilliard equation (also known as model B): 

d9(r,t) . 5F , . . , 

where e is a Langevin stochastic term which simulates thermal fluctuations as proposed by 



H.E. Cook in [p 



< e(r,t) > = (4.2) 
< e(r, t)e(r', t') > = -2k B TMA5(r - r')5(t - t') (4.3) 



In the previous Eq. |4.1| , M is a mobility constant which value may be related to the 
Fick diffusion constant Dpick- As known, for a given chemical species the parameter Dpick 
strongly depends on both temperature and coverage (see Ref. |HJ). For simplicity, we assum 
here that Dp ic k w 10~ 6 cm 2 / s (see ||16|| ). Neglecting the elastic term E e i, a linear expansion 
of Eq. O around a uniform coverage 8(r) = 0.5 shows that the Eq. [4.1| is equivalent 



to a diffusion equation for short wavelength perturbations (see ||12|| ). Then, we establish 
M.Fq.^ = D Fick which fixes approximately the mobility since ^£(#o — 0.5) = 16, in the 
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,2 ? 

present case. As both D Fick and |^ strongly decrease with a decreasing #0, the mobility 
variations with respect to the average coverage can not be estimated in the general context 
of our work. For simplicity, we choose to keep it constant, i.e., it is equal to the mobility 
estimated at half coverage. 

Regarding the kinetics of the coarse-grained surface, it is simply derived from the space 



and time discretization of Eq. |4.1| . The previous linear development of Eq. |4J] is still valid 
and it comes that the diffusion constant Dpick is rescaled as Dpick/d 2 due to the discretized 
Laplacian operator. It corresponds to the intuitive idea according to which the larger the 
coarse-graining is, i.e., the larger is d, the slower each unit cell of the coarse-graining evolves 
and thus the slower is the global kinetics. One guesses that the larger is a B domain, the 
longer is the time needed to transform its whole area in a A phase region. Thus we are able 
to propose a rough estimation of the surface kinetics time scale. Our computer calculations 
work with a time unit which is equivalent to 5t = d 2 /Dp ic k « 10 _8 s. To perform the time 
integration of Eq. |Q, we use a time increment of (10 -3 x 5t). 



The Eq. (3~T) is integrated starting from a uniform initial state for which 6{r) = 9q. 
Spontaneous formation of domains is observed for a coverage #0 ly m g between 0.25 < 6 < 
0.75. This range is imposed by the classical spinodal region defined by d 2 f / d 2 9 < (see [[0|). 
These limits are indeed slightly modified in presence of elastic interactions. The temperature 
dependency may be activated in the functional F c h em , in the diffusion coefficient Dp^, and 
in the amplitude of the Langevin noise. The external parameters are the material elastic 
constants that fix Xi the average coverage 80, the surface stress A and the a angle for both 
isotropic and anisotropic stress cases. 

Our numerical findings confirm the analytical results of Refs. |]^| : the phase separation 
is frustrated by long range interaction due to the crystal surface stress, i.e., the kinetics 
yields a final steady state. As model input parameters, one may retain that the domain 
boundary energy is of order / = lOmeV/A and the stress amplitude is either A = 40 mJ/m 2 
(or A = 0.25eV/A 2 ) for the Copper, Gold and Niobium crystal cases or A = 57 mJ/m 2 (or 
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A = 0.35eVyA 2 ) for the Chromium case. 

In the case of a broken surface symmetry with uniaxial stress, the kinetics yields non 
branched stripes structures with spacial periodicity along a given direction as shown on 
Fig.[V| (for a = 7r/4). The selected period strongly decreases when A increases but the 
formula proposed in Refs. which gives the period versus the model parameters, has not 
been confirmed by our computations because it required many long time computer runs. It 
is the purpose of a forthcoming paper. For a coverage 6q = 0.3 (see second row in Fig.[V|), 
the time needed to reach the final stripe structure is longer than for 8q = 0.5. As soon as the 
nucleation occures, it appears some distorted stripes length of which increases with coverage. 
During the coalescence regime, at 9q = 0.5, the stripes branch to each other before their 
connections move along the stripe sides such that it eliminates the domain which separates 
the primary stripe couple. The mecanism is very different at low coverage (see second row 
of Fig.[V|), since the short domains move to align themself with longer stripes. 

For a preserved surface symmetry, the kinetics computer calculations results are pre- 
sented on Figs.[V| and [V]. In Fig.[V], the time evolution of a half covered surface with and 
without induced stress is shown. One notes that for a negligible elastic interaction, i.e., 
A = 0, the domains never stop growing because of the Ostwald's ripening [0,0. We had 
to stop the simulation when the size of the domains is of the order of the whole sample. 
For A ^ 0, the kinetics is radically different since the system is driven to a "self-organized" 
final state. Beyong an evolving stage which duration decreases as A increases, the surface 
reaches a steady structure as shown by the two last pictures of the second and the third 
rows in Fig. |V|. The specific size and the distance between the neighboring domains are 
determined by the value of A. The larger is A, the smaller are the domain sizes. Accord- 
ing to our calculations, the final state is composed of branched stripes with tips as shown 
by the six last pictures in Fig.|V|. As in Fig.[V], the SO appears as soon as the nucleation 
occures. At this stage, the domains already arrange along preferential orientations. The 
structure evolves such that the thinest branches disappear to the benefit of the thick stripes 
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by first consuming the tips. The stripe orientations depend clearly on the sign of x- f° r 
X > as for Chromium, stripes are along either [110] or [110] while for x < as for Copper, 
stripes are along either [100] or [010]. According to the few experimental results about SO 
of chemisorbed mono-layer, the adsorption of Nitrogen on a (001) Chromium surface has 
been analysed with STM by M. Pinczolits et al. in Ref. ]T9|| . Infortunately, none of their 
results allows to conclude about the possible arrangement of Nitrogen atoms along [110] or 
[110] directions. The experimental results obtained by ||[7] about (001) Copper surface are 
in a good agreement with our calculations. 

In Fig.[V], the plot of the ratio SF/(Fq.S) where S is the total surface size shows how 
the relative free energy AF evolves with time. This quantity SF is the total free energy 
measured with respect to the free energy of a surface with only two compact domains of 
A and B phases. The stressed surface decomposition kinetics generates a steady state with 
lower energy because of a negative elastic energy E e \. After a nucleation regime where the 
total energy decreases strongly as time increases, the energy reaches an asymptotic value. 
The asymptotic state free energy clearly depends on the coverage 6 and the difference 
5F(8q = 0.5) — SF(9 = 0.25) increases in absolute value with A. The same plot with 
no surface stress shows that the surface evolves in order to reduce its positive domain 
boundary energy. In such a case, the asymptotic energy difference of two surfaces with 
different coverages, i.e., SF(9 Q = 0.5) — 6F(9o = 0.25) becomes negligible. 

At low coverage, as it was predicted in Ref. ||, an assembly of dots appear in the final 
stage. As for stripes, the arrangement of the dots also depends on the sign of x- The dots 
are arranged in an assembly of 1-dimensional rafts: the raft orientation is either along [110] 
or [110] for x > (first row of FigjV]), and either along [100] or [010] for x < (second 
row of Fig.[V|). In Fig.[V], the dots may be square shaped because of the relaxation of their 
elastic energy and the underlaying crystal symmetry. Comparison of the total free energy 
of a single A phase dot on a B surface, imposing different shapes to this dot, shows that 
the square shape is optimal, provided the area covered by the dot is larger than a critical 
size. For smaller sizes, the droplets are disk shaped. This critical size depends on both 

13 



the elastic constants and A. As a consequence, in the first stage of the growth, when the 
domains have small area, all of them are disk shaped. In the case of Copper, before the 
surface reaches its steady state, most of the islands undergo a shape tranformation passing 
from a disk shape to a square shape (see the first row of Fig.|V|). In such a case, the domain 
size is clearly larger than the critical shape transformation size. In the left picture of the first 
row of Fig.|V|, performed with Chromium elastic coefficients, it is the opposite case, since 
the disk shape domains never transform. In the right picture of the same row, increasing 
the surface coverage 8q yields larger domains. Some of them overpasses the critical size to 
becomes square with [110] and [110] oriented sides. 

For a coverage lying between 8 = 0.25 and half coverage, we find that the nanostructure 
of the final stage is composed of a mixing between short length stripes and dots (see Fig. |V|). 
So, we deduce that the surface SO shows a crossover with respect to the coverage 9 , i.e., 
going from droplet structure at low coverage to a stripe structure at half coverage, passing 
through mixed structures. 

The Fig. [V] shows the displacement field involved by the presence of the adsorbate. 
Nearby the surface, vortices are induced inside the bulk and extend over a depth of few 
ten nanometers. The maximum amplitude of the displacements is of few 0.1 A. With 
both a molecular dynamics study and a linear elasticity analysis of the specific N/Cu (001) 
case, those structures are also found by B. Croset et al. [17|] who first pointed them to us. 
Moreother they showed that the bulk displacements due to vortices contribute to X-ray 
diffraction and they deduce a surface force amplitude measurement: 2.2 10 -9 N.at~l (or 0.3 
eV/A 2 ). In our model, the A amplitude has same order of magnitude. 

For (001) Copper surface case, an induced shear is tested, i.e., a non zero /i = 5a^ 2 = da^i 
(defined above). Both the island shape and the arrangement directions are tilted with an 
angle that depends on the ratio /x/A (see Fig. [V]). 

The previous computations are performed with no thermal fluctuations. Putting on the 



Langevin noise in Eq. (JO), with an amplitude which is related to a temperature by Eq. 



73D, we performed the same simulations and obtained similar final states with the same 
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structures. This demonstrates the stability of both the kinetics and the final state. At half 
coverage and for the Copper elastic coefficients, the amplitude of the thermal noise may be 
increased up to the melting point of the Copper metal. Our method does not simulate the 
desorption process which could be activated at lower temperature than the crystal melting 
point. Moreother the white noise amplitude is calculated to satisfy the fluctuation dissipation 
theorem which is strictly valid around the thermodynamical equilibrium. Nevertheless the 
use of the Langevin noise allows to estimate qualitatively the surface structure ability to 
stand up to thermal fluctuations. Here the nanometric self-organised surface structure are 
proved to be very strong. 

Similar results have been obtained with elastic constants of Gold (x < ) and with 
Niobium elastic constants (x > 0) with a smaller amplitude A, indeed (see above for details). 



V. CONCLUSION AND PERSPECTIVES 

In the present paper, we investigated the patterns that are yielded from the (001) cubic 
crystal surface SO. If the surface square symmetry is preserved, the nanometric morpholo- 
gies are proved to depend on both the coverage and the elastic coefficients of the crystal. 
One-dimentional rafts of dots (with disk shape or square shape) occur at low coverage and 
branched stripes appear at half coverage. Rafts and stripes are aligned with either [100] or 
[010] for a negative anisotropy factor and with either [110] or [110] for a positive anisotropy 
factor. Our computations have been performed with different cubic metal elastic coeffi- 
cients, e.g., Gold, Copper, Chromium and Niobium. Obviously, it is possible to extend our 
calculations to the surfaces of either alloys or ionic crystal with cubic symmetry. In case of 
metal, this study should be extended to other surface symmetries which exibits interesting 
properties such as the (111) Gold surface e.g., the Gold Herringbone reconstruction (see, 
e.g., Refs. P,pO|), and as the (110) Copper surface, e.g., the Oxyde adsorption (see Ref. 



2TJ| ) . Then, in the framework of the spinodal theory for SO, only the elastic Green function 



should be modified. 
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In some specific cases, the crystal surface SO displays experimental features that are not 
included in the present model. The STM analysis of N/Cu(001) performed by Ellmer et al. |7j 
demonstrates that the surface morphologies are not equivalent for coverage 8q = (0.5+2) and 
for #o = (0.5 — a;), as one would expect from the spinodal decomposition theory we proposed. 
In addition, attractive interactions between Nitrogen islands seems to be dominant at very 
low coverage. We also note a remarkable result of T.W. Fishlock et al. |§ on Br/Cu(001) 
where STM analysis shows bromine islands organized in a chessboard at half coverage while 
we would expect branched stripes. We believe that our simple method can be improved 
in order to investigate those specific cases even if it is clear that a continuous approach is 
limited to the mesoscopic scale analysis, e.g., it seems hopeless to include certain atomic 
scale features as surface dislocations, or ad-atom missing row. 



In the present work, the chemisorbed layers thickness is neglected. As shown in Ref. [22 



for the Silicon quantum dots, this thickness plays a role in the arrangement of dots but for 
the moment no idea emerges to tackle this problem on the basis of our approach. 

For simplicity, we ignored the possible internal structures of each phases. If not it may 
involve both anisotropy and shear that are expected to modify the mesoscopic patterns. 
Moreother, some translational and orientational variants may be yielded from the different 
way the ad-atoms arrange on the underlying surface crystal lattice. Those variants can be 
treated in the present model by developing the free energy density with respect to the so 
called long range order parameters. This extension is well known within the phase field 



method for alloy physics fL0| . As a consequence of the coexistent variants, some anti-phase 
boundaries should appear between domains and play a role in the pattern growth. 

In summary, the model we present here is a promising start to describe and also to 
predict the nanostructures induced by self-organization of a crystal surface. 
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Figure captions 

Fig. 1 (A-F): Phase separation kinetics of a (001) cubic crystal surface (x — XAu ~ —0.5) 
with a stress directional anisotropy a = 7r/4 (defined in the text). Pictures occurence time 
is indicated at the low left hand corners: time unit is around 10 ns. Kinetics is started at 
initial time t — with a uniform surface state. The direction [010] is indicated. 

Fig. 2 (A-I): Phase separation kinetics on a (001) cubic crystal surface with negligible 
stress, i.e., A = (defined in the text) first row, and with A = 40 mJ/m 2 x — Xcu — — 1 , 
second row and A = 57 m.J/m 2 and x — Xcr = 1 third rows. 

Fig. 3: Time evolution of the total free energy SF for negligible stress, i.e., A = 0, and 
for a stressed surface, i.e., A = 40 m.J/m 2 for different coverages: 9 = 0.25 (dashed lines) 
and #o = 0.5 (solid lines). Time unit is 10 ns and the free energy is divided by the constant 
F .S (defined in the text). 

Fig. 4 (A-B) and 5 (A-B): First row: steady states of a (001) cubic crystal surface with 
A = 40 mJ/m 2 , x = Xcu = -1, #o = 0.21 (on the left) and 6 = 0.35 (on the right). 
Second row: steady state of a (001) cubic crystal surface with A = 57 mJ/m 2 , x = Xcr = 1 
, # = 0.21 (on the left) and 9 = 0.35 (on the right). 

Fig. 6 : Profile of the displacement field inside the bulk, for A = 40 mJ/m 2 , x = Xcu — 
— 1 and 9 = 0.21. The arrows lenght is proportional to the displacement (multiplied by 
one hundred). A distance d = lnm separates neighboring sites of the coarse-grained lattice. 
The directions [010] and [001] are indicated. 

Fig. 7: Steady state of a (001) cubic crystal surface with A = 40 mJ/m 2 , x = Xcu = — 1 
, #o = 0.3 and a shear component Sa^ 2 = <5ci2 = A/2. 
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